library(tidyverse)
library(RColorBrewer)
skyblue <- rgb(236, 240, 241, max = 255)

covered <- function(x, se, mu, q_z) {
  ci_hi <- x + q_z * se
  ci_lo <- x - q_z * se
  as.numeric(ci_hi > mu & ci_lo < mu)
}

col_spec <- cols(.default = col_double(), dgp = col_character())
fn <- list.files("output/", pattern = "^(?:main|dense|largen).*.csv")
out <- lapply(
  fn,
  function(x) read_csv(paste0("output/", x), col_types = col_spec)
)
out <- bind_rows(out)
out <- out %>%
  mutate(
    q_95 = qt(p = 0.975, n_units - n_covs),
    lasso_cov = covered(lasso_int, lasso_intse, 1, q_95),
    pds_cov = covered(pds_int, pds_intse, 1, q_95),
    oracle_cov = covered(oracle_int, oracle_intse, 1, q_95),
    fully_cov = covered(fully_int, fully_intse, 1, q_95),
    single_cov = covered(single_int, single_intse, 1, q_95)
  ) %>%
  pivot_longer(cols = oracle_main:single_cov,
               names_to = c("method", "qoi"),
               names_sep = "_")


## Absolute Bias
res <- out %>% filter(qoi == "int") %>%
  filter(!(method == "fully" & n_covs == 200 & n_units == 425)) %>% 
  group_by(dgp, n_units, n_covs, R2_y, R2_d, method) %>%
  summarize(bias = abs(mean(value - 1)),
            rmse = sqrt(mean((value - 1)^2)),
            se = sd(value),
            mean_est = mean(value))
res

res_se <- out %>%
  filter(qoi == "intse")  %>%
  filter(!(method == "fully" & n_covs == 200 & n_units == 425)) %>%
  group_by(dgp, n_units, n_covs, R2_y, R2_d, method) %>%
  summarize(seavg = mean(value))
res_se
res_se <- merge(res, res_se)
res_se$se_bias <-  res_se$seavg - res_se$se

res_cov <- out %>%
  filter(qoi == "cov") %>%
  filter(!(method == "fully" & n_covs == 200 & n_units == 425)) %>%  
  group_by(dgp, n_units, n_covs, R2_y, R2_d, method) %>%
  summarize(coverage = mean(value))
res_cov

## themes
theme_set(theme_minimal() +
            theme(panel.border = element_rect(fill = NA)))
theme_cousteau <- theme_light() +
  theme(plot.background = element_rect(fill = skyblue, color = skyblue),
        legend.background = element_rect(fill = skyblue))



meth_labs <- c(alasso = "ALasso (CV)",
               bart = "BART",
               fully = "Fully Moderated",
               krls = "KRLS",
               lasso = "Post-Lasso",
               oracle = "Oracle",
               pds = "Post-Double Selection",
               single = "Single Interaction")
hues <- seq(15, 375, length = 7 + 1)
cols <-  hcl(h = hues, l = 65, c = 100)[c(1, 2, 2, 3, 4, 5, 6, 7)]## viridisLite::viridis(n = 7)[c(1, 2, 2, 3, 4, 5, 5, 6, 6, 7)]
names(cols) <- c("pds", "lasso", "alasso", "oracle", "fully",
                 "krls", "bart", "single")
shapes <- c(15, 18, 18, 3, 17, 7, 19, 9)
names(shapes) <- c("pds", "lasso", "alasso", "oracle", "fully",
                 "krls", "bart", "single")

base_meths <- c("alasso", "fully", "krls", "oracle", "pds")
base_labs <- meth_labs[base_meths]
base_cols <- cols[base_meths]
base_shapes <- shapes[base_meths]

all_meths <- c("alasso", "fully", "bart", "krls", "pds", "single")
all_labs <- meth_labs[all_meths]
all_cols <- cols[all_meths]
all_shapes <- shapes[all_meths]

alt_meths <- c("fully", "lasso", "pds", "oracle")
alt_labs <- meth_labs[alt_meths]
alt_cols <- cols[alt_meths]
alt_shapes <- shapes[alt_meths]

facet_labs <- label_bquote(rows = .(n_covs)~"Covariates",
                           cols = {R^2}[y] == .(R2_y))
xlab <- expression(
  paste(
    "Partial ", R ^ 2, " for ", X[i] * V[i], " interaction on ", D[i]
  )
)

rmse_base <- res %>% filter(dgp == "main" & n_units == 425) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab,
       y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) 

rmse_base_all <- res %>%
  filter(dgp == "main" & n_units == 425) %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

rmse_base_alt <- res %>%
  filter(dgp == "main" & n_units == 425) %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base <- res %>%
  filter(dgp == "main" & n_units == 425) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method, shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_all <- res %>%
  filter(dgp == "main" & n_units == 425) %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_alt <- res %>%
  filter(dgp == "main" & n_units == 425) %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))


sebias_base <- res_se %>%
  filter(dgp == "main" & n_units == 425) %>%
  ggplot(aes(x = R2_d, y = se_bias, group = method, color = method,
             shape = method)) +
  labs(title = "Bias of SEs",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Bias",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = meth_labs, values = cols) +
  scale_shape_discrete(labels = meth_labs) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

cov_base <- res_cov %>%
  filter(dgp == "main" & n_units == 425) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)

cov_alt <- res_cov %>%
  filter(dgp == "main" & n_units == 425) %>%  
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)



## main RMSE
cairo_pdf("figures/main-rmse-sim.pdf", width = 8, height = 5, family = "Verdana")
rmse_base + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-rmse-sim-all.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-rmse-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-bias-sim.pdf", width = 8, height = 5, family = "Verdana")
bias_base + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-bias-sim-all.pdf", width = 8, height = 5, family = "Verdana")
bias_base_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-bias-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
bias_base_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-cov-sim.pdf", width = 8, height = 5, family = "Verdana")
cov_base + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/main-cov-alt-sim.pdf", width = 8, height = 5, family = "Verdana")
cov_alt + geom_point(size = 3) + geom_line()
dev.off()


## dense DGP
rmse_base_dense <- res %>% filter(dgp == "dense") %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) 

rmse_base_dense_all <- res %>%
  filter(dgp == "dense") %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

rmse_base_dense_alt <- res %>%
  filter(dgp == "dense") %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_dense <- res %>%
  filter(dgp == "dense") %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method, shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_dense_all <- res %>%
  filter(dgp == "dense") %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_dense_alt <- res %>%
  filter(dgp == "dense") %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))


sebias_base_dense <- res_se %>%
  filter(dgp == "dense") %>%
  ggplot(aes(x = R2_d, y = se_bias, group = method, color = method,
             shape = method)) +
  labs(title = "Bias of SEs",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Bias",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = meth_labs, values = cols) +
  scale_shape_discrete(labels = meth_labs) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

cov_base_dense <- res_cov %>%
  filter(dgp == "dense") %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)

cov_alt_dense <- res_cov %>%
  filter(dgp == "dense") %>%  
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)


cairo_pdf("figures/dense-rmse-sim.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_dense + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-rmse-sim-all.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_dense_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-rmse-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_dense_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-bias-sim.pdf", width = 8, height = 5, family = "Verdana")
bias_base_dense + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-bias-sim-all.pdf", width = 8, height = 5, family = "Verdana")
bias_base_dense_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-bias-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
bias_base_dense_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-cov-sim.pdf", width = 8, height = 5, family = "Verdana")
cov_base_dense + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/dense-cov-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
cov_alt_dense + geom_point(size = 3) + geom_line()
dev.off()



## large_n DGP
rmse_base_large_n <- res %>% filter(dgp == "main" & n_units == 1000) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) 

rmse_base_large_n_all <- res %>%
  filter(dgp == "main" & n_units == 1000) %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

rmse_base_large_n_alt <- res %>%
  filter(dgp == "main" & n_units == 1000) %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = rmse, group = method, color = method,
             shape = method)) +
  labs(title = "Root Mean Square Error",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "RMSE",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_large_n <- res %>%
  filter(dgp == "main" & n_units == 1000) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method, shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_large_n_all <- res %>%
  filter(dgp == "main" & n_units == 1000) %>%
  filter(method %in% all_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = all_labs, values = all_cols) +
  scale_shape_manual(labels = all_labs, values = all_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

bias_base_large_n_alt <- res %>%
  filter(dgp == "main" & n_units == 1000) %>%
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = bias, group = method, color = method,
               shape = method)) +
  labs(title = "Absolute Bias",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "|Bias|",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))


sebias_base_large_n <- res_se %>%
  filter(dgp == "main" & n_units == 1000) %>%
  ggplot(aes(x = R2_d, y = se_bias, group = method, color = method,
             shape = method)) +
  labs(title = "Bias of SEs",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Bias",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = meth_labs, values = cols) +
  scale_shape_discrete(labels = meth_labs) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

cov_base_large_n <- res_cov %>%
  filter(dgp == "main" & n_units == 1000) %>%  
  filter(method %in% base_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = base_labs, values = base_cols) +
  scale_shape_manual(labels = base_labs, values = base_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)

cov_alt_large_n <- res_cov %>%
  filter(dgp == "main" & n_units == 1000) %>%  
  filter(method %in% alt_meths) %>%
  ggplot(aes(x = R2_d, y = coverage, group = method, color = method,
             shape = method)) +
  labs(title = "Coverage of 95% Confidence Intervals",
       subtitle = "Quadratic Decay of Covariate Effects",
       x = xlab, y = "Coverage",
       color = "Method", shape = "Method") +
  facet_grid(n_covs ~ R2_y, labeller = facet_labs) +
  scale_color_manual(labels = alt_labs, values = alt_cols) +
  scale_shape_manual(labels = alt_labs, values = alt_shapes) +
  theme(axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) +
  geom_abline(intercept = 0.95, slope = 0, linetype = 2)


cairo_pdf("figures/large_n-rmse-sim.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_large_n + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-rmse-sim-all.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_large_n_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-rmse-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
rmse_base_large_n_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-bias-sim.pdf", width = 8, height = 5, family = "Verdana")
bias_base_large_n + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-bias-sim-all.pdf", width = 8, height = 5, family = "Verdana")
bias_base_large_n_all + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-bias-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
bias_base_large_n_alt + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-cov-sim.pdf", width = 8, height = 5, family = "Verdana")
cov_base_large_n + geom_point(size = 3) + geom_line()
dev.off()
cairo_pdf("figures/large_n-cov-sim-alt.pdf", width = 8, height = 5, family = "Verdana")
cov_alt_large_n + geom_point(size = 3) + geom_line()
dev.off()

